#pragma once

#include <cmath>
#include <vector>
#include <numeric>
#include <algorithm>

// 检测SIMD支持
#if defined(__AVX__)
    #include <immintrin.h>
    #define HFT_SIMD_AVX 1
#elif defined(__SSE4_1__)
    #include <smmintrin.h>
    #define HFT_SIMD_SSE4_1 1
#elif defined(__SSE3__)
    #include <pmmintrin.h>
    #define HFT_SIMD_SSE3 1
#elif defined(__SSE2__) || (defined(_MSC_VER) && (defined(_M_X64) || (defined(_M_IX86_FP) && _M_IX86_FP >= 2)))
    #include <emmintrin.h>
    #define HFT_SIMD_SSE2 1
#endif

namespace utils {
namespace simd {

/**
 * @brief 检查当前CPU是否支持SIMD指令集
 * @return 支持的SIMD级别
 */
enum class SimdLevel {
    None,
    SSE2,
    SSE3,
    SSE4_1,
    AVX
};

inline SimdLevel getSimdLevel() {
#if defined(HFT_SIMD_AVX)
    return SimdLevel::AVX;
#elif defined(HFT_SIMD_SSE4_1)
    return SimdLevel::SSE4_1;
#elif defined(HFT_SIMD_SSE3)
    return SimdLevel::SSE3;
#elif defined(HFT_SIMD_SSE2)
    return SimdLevel::SSE2;
#else
    return SimdLevel::None;
#endif
}

/**
 * @brief 使用SIMD指令计算向量的和
 * @param data 数据指针
 * @param size 数据大小
 * @return 向量元素之和
 */
inline double sum(const double* data, size_t size) {
#if defined(HFT_SIMD_AVX)
    // 使用AVX指令集
    const size_t simd_size = 4; // __m256d处理4个double
    size_t simd_count = size / simd_size;
    size_t remainder = size % simd_size;
    
    __m256d sum_vec = _mm256_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m256d data_vec = _mm256_loadu_pd(data + i * simd_size);
        sum_vec = _mm256_add_pd(sum_vec, data_vec);
    }
    
    // 水平求和
    __m128d sum_high = _mm256_extractf128_pd(sum_vec, 1);
    __m128d sum_low = _mm256_castpd256_pd128(sum_vec);
    __m128d sum_hl = _mm_add_pd(sum_high, sum_low);
    __m128d sum_lh = _mm_permute_pd(sum_hl, 1);
    __m128d sum_scalar = _mm_add_sd(sum_hl, sum_lh);
    
    double result = _mm_cvtsd_f64(sum_scalar);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        result += data[i];
    }
    
    return result;
#elif defined(HFT_SIMD_SSE2)
    // 使用SSE2指令集
    const size_t simd_size = 2; // __m128d处理2个double
    size_t simd_count = size / simd_size;
    size_t remainder = size % simd_size;
    
    __m128d sum_vec = _mm_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m128d data_vec = _mm_loadu_pd(data + i * simd_size);
        sum_vec = _mm_add_pd(sum_vec, data_vec);
    }
    
    // 水平求和
    __m128d sum_high = _mm_unpackhi_pd(sum_vec, sum_vec);
    __m128d result_vec = _mm_add_sd(sum_vec, sum_high);
    
    double result = _mm_cvtsd_f64(result_vec);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        result += data[i];
    }
    
    return result;
#else
    // 标量实现
    return std::accumulate(data, data + size, 0.0);
#endif
}

/**
 * @brief 使用SIMD指令计算向量的平方和
 * @param data 数据指针
 * @param size 数据大小
 * @return 向量元素的平方和
 */
inline double sumOfSquares(const double* data, size_t size) {
#if defined(HFT_SIMD_AVX)
    // 使用AVX指令集
    const size_t simd_size = 4; // __m256d处理4个double
    size_t simd_count = size / simd_size;
    size_t remainder = size % simd_size;
    
    __m256d sum_vec = _mm256_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m256d data_vec = _mm256_loadu_pd(data + i * simd_size);
        sum_vec = _mm256_add_pd(sum_vec, _mm256_mul_pd(data_vec, data_vec));
    }
    
    // 水平求和
    __m128d sum_high = _mm256_extractf128_pd(sum_vec, 1);
    __m128d sum_low = _mm256_castpd256_pd128(sum_vec);
    __m128d sum_hl = _mm_add_pd(sum_high, sum_low);
    __m128d sum_lh = _mm_permute_pd(sum_hl, 1);
    __m128d sum_scalar = _mm_add_sd(sum_hl, sum_lh);
    
    double result = _mm_cvtsd_f64(sum_scalar);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        result += data[i] * data[i];
    }
    
    return result;
#elif defined(HFT_SIMD_SSE2)
    // 使用SSE2指令集
    const size_t simd_size = 2; // __m128d处理2个double
    size_t simd_count = size / simd_size;
    size_t remainder = size % simd_size;
    
    __m128d sum_vec = _mm_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m128d data_vec = _mm_loadu_pd(data + i * simd_size);
        sum_vec = _mm_add_pd(sum_vec, _mm_mul_pd(data_vec, data_vec));
    }
    
    // 水平求和
    __m128d sum_high = _mm_unpackhi_pd(sum_vec, sum_vec);
    __m128d result_vec = _mm_add_sd(sum_vec, sum_high);
    
    double result = _mm_cvtsd_f64(result_vec);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        result += data[i] * data[i];
    }
    
    return result;
#else
    // 标量实现
    return std::inner_product(data, data + size, data, 0.0);
#endif
}

/**
 * @brief 使用SIMD指令计算向量的均值
 * @param data 数据指针
 * @param size 数据大小
 * @return 向量元素的均值
 */
inline double mean(const double* data, size_t size) {
    if (size == 0) return 0.0;
    return sum(data, size) / size;
}

/**
 * @brief 使用SIMD指令计算向量的标准差
 * @param data 数据指针
 * @param size 数据大小
 * @return 向量元素的标准差
 */
inline double stdDev(const double* data, size_t size) {
    if (size <= 1) return 0.0;
    
    double mean_val = mean(data, size);
    double sum_sq = 0.0;
    
#if defined(HFT_SIMD_AVX)
    // 使用AVX指令集
    const size_t simd_size = 4; // __m256d处理4个double
    size_t simd_count = size / simd_size;
    
    __m256d mean_vec = _mm256_set1_pd(mean_val);
    __m256d sum_vec = _mm256_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m256d data_vec = _mm256_loadu_pd(data + i * simd_size);
        __m256d diff_vec = _mm256_sub_pd(data_vec, mean_vec);
        sum_vec = _mm256_add_pd(sum_vec, _mm256_mul_pd(diff_vec, diff_vec));
    }
    
    // 水平求和
    __m128d sum_high = _mm256_extractf128_pd(sum_vec, 1);
    __m128d sum_low = _mm256_castpd256_pd128(sum_vec);
    __m128d sum_hl = _mm_add_pd(sum_high, sum_low);
    __m128d sum_lh = _mm_permute_pd(sum_hl, 1);
    __m128d sum_scalar = _mm_add_sd(sum_hl, sum_lh);
    
    sum_sq = _mm_cvtsd_f64(sum_scalar);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        double diff = data[i] - mean_val;
        sum_sq += diff * diff;
    }
#elif defined(HFT_SIMD_SSE2)
    // 使用SSE2指令集
    const size_t simd_size = 2; // __m128d处理2个double
    size_t simd_count = size / simd_size;
    
    __m128d mean_vec = _mm_set1_pd(mean_val);
    __m128d sum_vec = _mm_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m128d data_vec = _mm_loadu_pd(data + i * simd_size);
        __m128d diff_vec = _mm_sub_pd(data_vec, mean_vec);
        sum_vec = _mm_add_pd(sum_vec, _mm_mul_pd(diff_vec, diff_vec));
    }
    
    // 水平求和
    __m128d sum_high = _mm_unpackhi_pd(sum_vec, sum_vec);
    __m128d result_vec = _mm_add_sd(sum_vec, sum_high);
    
    sum_sq = _mm_cvtsd_f64(result_vec);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        double diff = data[i] - mean_val;
        sum_sq += diff * diff;
    }
#else
    // 标量实现
    for (size_t i = 0; i < size; ++i) {
        double diff = data[i] - mean_val;
        sum_sq += diff * diff;
    }
#endif
    
    return std::sqrt(sum_sq / (size - 1));
}

/**
 * @brief 使用SIMD指令计算两个向量的相关系数
 * @param data1 第一个向量的数据指针
 * @param data2 第二个向量的数据指针
 * @param size 数据大小
 * @return 相关系数
 */
inline double correlation(const double* data1, const double* data2, size_t size) {
    if (size <= 1) return 0.0;
    
    double mean1 = mean(data1, size);
    double mean2 = mean(data2, size);
    
    double sum_xy = 0.0;
    double sum_x2 = 0.0;
    double sum_y2 = 0.0;
    
#if defined(HFT_SIMD_AVX)
    // 使用AVX指令集
    const size_t simd_size = 4; // __m256d处理4个double
    size_t simd_count = size / simd_size;
    
    __m256d mean1_vec = _mm256_set1_pd(mean1);
    __m256d mean2_vec = _mm256_set1_pd(mean2);
    __m256d sum_xy_vec = _mm256_setzero_pd();
    __m256d sum_x2_vec = _mm256_setzero_pd();
    __m256d sum_y2_vec = _mm256_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m256d data1_vec = _mm256_loadu_pd(data1 + i * simd_size);
        __m256d data2_vec = _mm256_loadu_pd(data2 + i * simd_size);
        
        __m256d diff1_vec = _mm256_sub_pd(data1_vec, mean1_vec);
        __m256d diff2_vec = _mm256_sub_pd(data2_vec, mean2_vec);
        
        sum_xy_vec = _mm256_add_pd(sum_xy_vec, _mm256_mul_pd(diff1_vec, diff2_vec));
        sum_x2_vec = _mm256_add_pd(sum_x2_vec, _mm256_mul_pd(diff1_vec, diff1_vec));
        sum_y2_vec = _mm256_add_pd(sum_y2_vec, _mm256_mul_pd(diff2_vec, diff2_vec));
    }
    
    // 水平求和
    __m128d sum_xy_high = _mm256_extractf128_pd(sum_xy_vec, 1);
    __m128d sum_xy_low = _mm256_castpd256_pd128(sum_xy_vec);
    __m128d sum_xy_hl = _mm_add_pd(sum_xy_high, sum_xy_low);
    __m128d sum_xy_lh = _mm_permute_pd(sum_xy_hl, 1);
    __m128d sum_xy_scalar = _mm_add_sd(sum_xy_hl, sum_xy_lh);
    
    __m128d sum_x2_high = _mm256_extractf128_pd(sum_x2_vec, 1);
    __m128d sum_x2_low = _mm256_castpd256_pd128(sum_x2_vec);
    __m128d sum_x2_hl = _mm_add_pd(sum_x2_high, sum_x2_low);
    __m128d sum_x2_lh = _mm_permute_pd(sum_x2_hl, 1);
    __m128d sum_x2_scalar = _mm_add_sd(sum_x2_hl, sum_x2_lh);
    
    __m128d sum_y2_high = _mm256_extractf128_pd(sum_y2_vec, 1);
    __m128d sum_y2_low = _mm256_castpd256_pd128(sum_y2_vec);
    __m128d sum_y2_hl = _mm_add_pd(sum_y2_high, sum_y2_low);
    __m128d sum_y2_lh = _mm_permute_pd(sum_y2_hl, 1);
    __m128d sum_y2_scalar = _mm_add_sd(sum_y2_hl, sum_y2_lh);
    
    sum_xy = _mm_cvtsd_f64(sum_xy_scalar);
    sum_x2 = _mm_cvtsd_f64(sum_x2_scalar);
    sum_y2 = _mm_cvtsd_f64(sum_y2_scalar);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        double diff1 = data1[i] - mean1;
        double diff2 = data2[i] - mean2;
        sum_xy += diff1 * diff2;
        sum_x2 += diff1 * diff1;
        sum_y2 += diff2 * diff2;
    }
#elif defined(HFT_SIMD_SSE2)
    // 使用SSE2指令集
    const size_t simd_size = 2; // __m128d处理2个double
    size_t simd_count = size / simd_size;
    
    __m128d mean1_vec = _mm_set1_pd(mean1);
    __m128d mean2_vec = _mm_set1_pd(mean2);
    __m128d sum_xy_vec = _mm_setzero_pd();
    __m128d sum_x2_vec = _mm_setzero_pd();
    __m128d sum_y2_vec = _mm_setzero_pd();
    
    // 使用SIMD处理大部分数据
    for (size_t i = 0; i < simd_count; ++i) {
        __m128d data1_vec = _mm_loadu_pd(data1 + i * simd_size);
        __m128d data2_vec = _mm_loadu_pd(data2 + i * simd_size);
        
        __m128d diff1_vec = _mm_sub_pd(data1_vec, mean1_vec);
        __m128d diff2_vec = _mm_sub_pd(data2_vec, mean2_vec);
        
        sum_xy_vec = _mm_add_pd(sum_xy_vec, _mm_mul_pd(diff1_vec, diff2_vec));
        sum_x2_vec = _mm_add_pd(sum_x2_vec, _mm_mul_pd(diff1_vec, diff1_vec));
        sum_y2_vec = _mm_add_pd(sum_y2_vec, _mm_mul_pd(diff2_vec, diff2_vec));
    }
    
    // 水平求和
    __m128d sum_xy_high = _mm_unpackhi_pd(sum_xy_vec, sum_xy_vec);
    __m128d sum_xy_result = _mm_add_sd(sum_xy_vec, sum_xy_high);
    
    __m128d sum_x2_high = _mm_unpackhi_pd(sum_x2_vec, sum_x2_vec);
    __m128d sum_x2_result = _mm_add_sd(sum_x2_vec, sum_x2_high);
    
    __m128d sum_y2_high = _mm_unpackhi_pd(sum_y2_vec, sum_y2_vec);
    __m128d sum_y2_result = _mm_add_sd(sum_y2_vec, sum_y2_high);
    
    sum_xy = _mm_cvtsd_f64(sum_xy_result);
    sum_x2 = _mm_cvtsd_f64(sum_x2_result);
    sum_y2 = _mm_cvtsd_f64(sum_y2_result);
    
    // 处理剩余元素
    for (size_t i = simd_count * simd_size; i < size; ++i) {
        double diff1 = data1[i] - mean1;
        double diff2 = data2[i] - mean2;
        sum_xy += diff1 * diff2;
        sum_x2 += diff1 * diff1;
        sum_y2 += diff2 * diff2;
    }
#else
    // 标量实现
    for (size_t i = 0; i < size; ++i) {
        double diff1 = data1[i] - mean1;
        double diff2 = data2[i] - mean2;
        sum_xy += diff1 * diff2;
        sum_x2 += diff1 * diff1;
        sum_y2 += diff2 * diff2;
    }
#endif
    
    double denominator = std::sqrt(sum_x2 * sum_y2);
    if (denominator < 1e-10) return 0.0;
    
    return sum_xy / denominator;
}